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First  and  Second  Year  Report  Period  (Dec.  2005  to  Dec.  2007) 


This  report  combines  the  first  two  years  of  the  formal  grant  award  period.  Because  of  a  Faculty 
Development  Leave  for  the  PI  and  the  associated  extreme  difficulty  in  identifying  and  hiring 
personnel  to  work  on  the  project,  it  was  impossible  to  achieve  the  planned  progress  for  the  first 
year.  This  difficulty,  as  well  as  the  inability  to  use  the  funds  for  the  planned  instrumentation 
(since  a  commercial  package  could  not  be  obtained  ),  led  Dr.  Goldwasser  and  the  PI  to  conclude 
that  a  no-cost  extension  of  the  grant  would  be  both  inevitable  and  necessary. 

2006 

The  PI  left  for  his  scheduled  year  of  Faculty  Developmental  Leave  in  January,  2006  beginning 
with  five  months  as  a  visiting  professor  at  the  University  of  Milan.  The  PI  returned  in  June, 
2006.  The  last  five  months  of  the  year  was  spent  as  a  Senior  Fulbright  Fellow  at  the  Wroclaw 
University  of  Technology.  Because  the  departure  for  Milan  was,  for  practical  purposes,  at  the 
inception  of  the  grant,  the  search  for  postdoctoral  associates  had  to  be  initiated  from  abroad.  A 
great  deal  of  time  was  taken  in  trying  to  locate  suitable  postdoctoral  associates  over  that  time. 
However,  few  with  the  necessary  background  were  identified  and  those  that  were,  either  had 
prior  commitments,  or  we  not  interested  in  pursuing  the  position  due,  in  no  small  part,  to  visa 
problems.  This  proved  to  be  quite  time-consuming  even  though  unsuccessful.  However,  with 
the  graduation  of  one  of  the  Pi’s  graduate  students,  who  had  worked  on  energetic  materials  for 
his  PhD,  it  appeared  that  it  could  be  possible  to  make  some  progress  on  the  research.  The 
individual,  J.  J.  Haycraft,  agreed  to  work  on  the  project.  He  began  to  devise  an  experimental 
protocol  for  measurement  and  introduction  of  defects  in  the  energetic  materials.  The  initial  work 
was  focused  on  design  and  construction  of  an  apparatus  for  injecting  defects  into  the  crystals 
using  PZT  ceramics.  Unfortunately,  Dr.  Haycraft  was  made  an  extremely  good  job  offer  arising 
from  a  prior  application.  He  left  in  August  for  his  new  position. 

During  this  time,  the  PI  was  engaged  in  working  with  his  collaborator,  Prof.  Gavezzotti  at  the 
University  of  Milan,  on  trying  to  develop  a  calculational  approach  to  understanding  the  role  of 
defects  in  solids  of  energetic  materials.  In  the  course  of  this  work,  it  was  found  that  the  density 
of  crystals  of  energetic  materials  are  abnormally  high  compared  to  crystals  of  other  organic 
molecular  crystals.  This  was  investigated  by  calculation  and  it  was  found  that,  contrary  to 
commonly  held  belief,  the  intermolecular  forces  in  these  solids  were  significantly  more  strong 
than  the  Coulombic  interactions.  The  results  of  these  lattice  energy  investigations  were 
submitted  to  the  J.  Phys.  Chem.  and  subsequently  published  in  that  Journal.  (  C.  J.  Eckhardt,  A. 
Gavezzotti,  J.  Phys.  Chem.  Bill  (2007)  3430.) 


2007 

With  Dr.  Haycraft’s  departure,  the  search  for  the  second  postdoctoral  associate  had  to  be 
expanded  to  a  search  for  two  associates.  Dr.  Himansu  Mohapatra  was  hired  to  continue  with  the 
work  begun  by  Dr.  Haycraft.  In  addition,  Dr.  Mohapatra  was  charged  with  the  growth  of 
suitable  crystals  of  the  less  common  polymorphs  of  HMX  and  RDX.  While  they  could  be 
obtained,  they  were  not  of  the  optical  quality  needed  for  the  proposed  experiments.  Dr. 
Mohapatra  spent  the  bulk  of  his  time  in  trying  to  devise  suitable  methods  of  growth  of  the 


polymorphs  of  the  two  energetic  materials.  The  search  continued  for  the  second  postdoctoral 
associate.  The  search  focused  more  on  scientists  in  the  US  in  order  to  avoid  the  problems  with 
foreign  workers.  Further,  it  became  evident  that  the  criteria  for  the  positions  would  have  to  be 
severely  modified.  This  permitted  identification  of  several  potential  people.  One,  in  particular, 
was  almost  perfectly  suited  to  the  position  but  declined  the  offer  in  favor  of  another.  This 
pattern  repeated.  However,  Dr.  Haidong  Zhang,  with  a  new  PhD  from  the  University  of  Florida, 
agreed  to  join  the  group  in  August,  2007.  This  was  particularly  welcome  since  Dr.  Zhang  had 
some  familiarity  with  working  at  synchrotrons.  The  light  source  for  the  proposed  vacuum 
ultraviolet  (VUV)  spectroreflectometer  is  to  be  the  CAMD  synchrotron  light  source  in  Baton 
Rouge. 

Upon  arrival,  Dr.  Zhang  began  immediately  to  work  on  the  design  of  the  VUV 
spectroreflectometer.  This  is  not  a  simple  task  since  the  sample,  an  organic,  is  inherently 
unstable  in  a  vacuum.  Nevertheless,  as  of  December,  2007,  Dr.  Zhang  had  made  very  good 
progress  in  identifying  the  problems  inherent  with  this  difficult  experiment.  By  December, 

2007,  he  was  identifying  what  components  were  viable  and  had  design  drawings  for  the 
instrument.  It  is  anticipated  it  will  be  installed  in  the  summer  of  2008.  Dr.  Zhang’s 
appointment  runs  to  August,  2008  and  he  has  agreed  to  stay  another  year  to  complete  the 
research. 

Near  the  end  of  2007,  a  second  postdoctoral  associate  with  experience  in  Brillouin  scattering  was 
identified.  Dr.  Zhando  Utegulov,  formerly  at  NIST,  agreed  to  join  the  project  and  work  on  the 
defect  problem.  Dr.  Mohapatra  had  indicated  that  he  would  rather  have  another  project  and  he 
has  been  able  to  successfully  help  Dr.  Utegulov  move  into  the  defect  studies  of  the  proposed 
research.  It  is  expected  that  the  proposed  studies  will  be  underway  by  the  summer  of  2008. 

Because  of  the  Pi’s  Faculty  Development  Leave  and  the  extreme  difficulty  in  staffing  for  this 
project,  it  has  been  effectively  delayed  by  at  least  a  year  and  perhaps  somewhat  more.  Now,  with 
full  staffing,  we  are  desperately  trying  to  make  up  for  lost  time.  The  project  is  now  well 
underway  and  significant  results  are  expected  by  the  end  of  2008. 

The  PI  applied  for,  and  was  granted,  the  Faculty  Leave  of  Absence  before  the  particular 

program  for  this  ONR  research  was  announced. 
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A  database  of  43  literature  X-ray  crystal  structure  determinations  for  compounds  with  known,  or  possible, 
energetic  properties  has  been  collected  along  with  some  sublimation  enthalpies.  A  statistical  study  of  these 
crystal  structures,  when  compared  to  a  sample  of  general  organic  crystals,  reveals  a  population  of  anomalously 
short  intermolecular  oxygen-oxygen  separations  with  an  average  crystal  packing  coefficient  of  0.77  that 
differs  significantly  from  0.70  found  for  the  general  population.  For  the  calculation  of  lattice  energies,  three 
atom— atom  potential  energy  schemes  and  the  semiempirical  SCDS-PIXEL  scheme  are  compared.  The  nature 
of  the  packing  forces  in  these  energetic  materials  is  further  analyzed  by  a  study  of  the  dispersive  versus 
Coulombic  contributions  to  overall  lattice  energies  and  to  molecule— molecule  energies  in  pairs  of  near 
neighbors  in  the  crystals,  a  partitioning  made  possible  by  the  unique  features  of  the  SCDS-PIXEL  scheme. 
It  is  shown  that  dispersion  forces  are  stronger  than  Coulombic  forces,  contrary  to  common  belief.  The  low 
abundance  of  hydrogen  atoms  in  these  molecules,  the  close  oxygen— oxygen  contacts,  and  the  high  packing 
coefficients  explain  the  observation  that,  for  these  energetic  materials,  crystal  densities  are  anomalously  high 
compared  to  those  of  most  organic  materials.  However,  an  understanding,  not  to  mention  prediction  or  control, 
of  the  deeper  mechanisms  for  the  explosive  power  of  these  crystalline  materials,  such  as  the  role  of  lattice 
defects,  remains  beyond  present  capabilities. 


I.  Introduction 

Weighing  the  advantages  and  dangers  of  developing  more 
efficient  energetic  materials  is  a  difficult  matter  that  involves 
ethics  and  politics  rather  than  pure  science,  but  it  requires  due 
consideration  of  the  fact  that  the  preponderant  use  of  explosives 
is  for  non-military  applications.  There  are  two  overriding 
concerns  for  those  who  use  these  materials:  safety  and  explosive 
power.  Energetic  materials  that  are  resistant  to  detonation  can 
be  handled  more  safely,  although  often  at  the  expense  of 
explosive  power.  These  matters  are  amenable  to  the  usual  studies 
of  structure— property  relationships,  hence  the  present  work. 

For  decades,  explosive  power  has  been  correlated  with  the 
energy  density  of  the  material  and,  in  particular,  with  the  mass 
density  because  of  the  dependence  of  the  magnitude  of  the 
detonation  velocity  and  of  the  detonation  sensitivity  to  this 
simple  physical  property.1  For  crystals  comprised  only  of  atoms 
of  elements  of  the  first  and  second  periods,  the  densities  of 
energetic  materials  tend  to  be  anomalously  large.  This  may  have 
implications  regarding  defects,  since  they  usually  lower  the 
density.  Another  important  measure  of  explosive  power  has  been 
the  oxygen  balance  that  is  taken  to  be  optimal  when  the  amount 
of  oxygen  in  the  energetic  material  is  just  sufficient  for 
conversion  of  all  oxygen  to  only  carbon  dioxide  and  water.2 

Most  energetic  materials  are  organic  molecules  of  rather 
irregular  shape,  and  thus  a  question  of  interest  is  how  these 
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molecules  pack  in  crystals  to  obtain  the  most  efficient  space 
occupation  and  hence  a  high  density  for  a  given  chemical 
composition.  For  such  oxygen-rich  substances,  it  would  nor¬ 
mally  be  expected  that  closely  packed  molecules  that  have  many 
peripheral  oxygen  atoms  would  find  difficulty  in  avoiding  close 
and,  presumably,  repulsive  contacts.  A  crucial  matter  then  is 
the  nature  and  magnitude  of  the  intermolecular  forces  that  bind 
the  molecules  in  these  crystals.  A  related  topic  is  crystal 
polymorphism,  because  different  polymorphs  exhibit  varying 
sensitivity  to  detonation,  thus  confirming  that  solid-state  proper¬ 
ties  must  have  significant  influence  on  the  mechanism  of 
detonation.18  A  proper  understanding  of  crystal  packing  forces 
may  help  in  controlling  the  growth  of  favorable  polymorphs. 
Further,  the  study  of  “perfect”  crystals  is  a  first  step  in  the 
analysis  of  the  defects  which  disrupt  the  crystal  lattice  and 
presumably  have  a  significant  influence  on  the  performance  of 
these  materials.3 

For  these  reasons,  we  have  undertaken  a  detailed  packing 
analysis  and  examination  of  the  binding  forces  of  a  series  of 
crystals  of  energetic  materials  by  a  combination  of  recently 
developed  statistical  techniques  and  molecular  simulation 
methods.4  Using  crystallographic  databases,  some  geometrical 
features  of  the  crystal  packing  of  energetic  materials  are 
compared  with  those  of  a  general  collection  of  organic 
compounds.  For  the  calculation  of  intermolecular  energies,  the 
available  empirical  atom— atom  force  fields  are  critically 
examined,  but  the  results  are  also  compared  with  those  obtained 
by  more  accurate  semiempirical  and  ab  initio  methods.  The 
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prospects  for  an  accurate  and  reliable  simulation  of  these 
materials  are  then  examined. 

II.  Experimental  Data 

1.  Heats  of  Sublimation.  The  sublimation  enthalpy,  A//sub, 
is  a  vital  piece  of  information  in  crystal -structure  analysis  and 
force-field  calibrations,  because  its  negative  should  be  equal  to 
the  calculated  lattice  energy.  The  NIST  thermodynamic  data¬ 
base5  was  searched  for  thermodynamic  properties  of  the 
crystalline  materials  of  interest  here,  in  particular  for  sublimation 
enthalpies  of  compounds  whose  crystal  structure  also  has  been 
determined  by  X-ray  diffraction.  The  primary  literature6  was 
also  searched  independently,  and  often  checked  to  help  in 
making  a  choice  between  different  A//SUb  values  stored  in  the 
NIST  database.  In  particular,  the  A7/SUb  for  MMX  ( 1 ,3,5,7- 
tetranitro-l,3,5,7-tetraazacyclooctane,  refcode  OCHTET)  has 
been  assessed  from  a  Clausius-Clapeyron  treatment  of  the 
original  data.60  The  reliability  of  thermochemical  data  is 
somewhat  variable,  and  discrepancies  of  up  to  20-30%  among 
different  values  for  A//sub  of  the  same  crystal  are  quite  common 
(for  an  extreme  instance,  the  reported  A//SUb’s  for  dimethyl- 
nitramine  range  from  41  to  70  kJ  mol1).  The  experimental  data 
come  from  different  techniques,  or  often  refer  to  high  temper¬ 
ature,  as  is  appropriate  for  the  measurement  of  vapor  pressures 
of  nonvolatile  compounds,  while  calculated  lattice  energies 
derive  from  structures  determined  at  room  temperature.  As¬ 
suming  a  typical  range  of  values  for  gas-crystal  ACp’s  for 
organic  crystals  of  20— 40  J  K_I  mol1,  one  can  estimate  a  2— 4 
kJ  mol-1  increase  in  the  A//SUb  for  a  100  K  lowering  of  the 
temperature.  In  view  of  the  inherent  uncertainties,  the  compari¬ 
son  between  experimental  A//SUb’s  and  calculated  lattice  energies 
is  no  more  than  a  guideline  in  the  calibration  of  crystal 
potentials,  and  only  an  overall  consideration  of  trends  over  many 
data  can  be  significant;  parametrization  based  on  only  a  few 
experimental  data  can  sometimes  be  misleading. 

2.  Crystal  Structures.  The  Cambridge  Structural  Database7 
(CSD)  and  the  primary  literature  were  searched  for  X-ray  crystal 
structure  determinations  of  materials  with  promise  for  high- 
energy  content,  especially  molecules  containing  several  nitro 
groups.  The  final  database  of  43  selected  crystals  (see  Tables  1 
and  2)  is  representative  rather  than  exhaustive,  but  proved  to 
be  sufficient  for  the  present  purposes.  The  crystal  structure 
determinations  are  of  variable  accuracy  as  a  consequence  of 
changeable  crystal  quality  and  the  other  usual  experimental 
limitations  of  the  X-ray  work.  This  hardly  affects  the  atom- 
atom  force-field  calculations,  whose  potentials  are  isotropic  and 
radial,  and  rather  insensitive  to  minor  structural  variations. 
However,  even  a  minor  inaccuracy  in  atomic  positions  signifi¬ 
cantly  affects  the  ab  initio  quantum  chemical  calculations  or 
the  PIXEL  calculations  (see  below)  that  take  into  account  all 
of  the  electron  density  distribution  and  are  much  more  sensitive 
to  structural  details.  For  example,  unrealistic  intramolecular 
energy  differences  are  calculated  ab  initio  for  the  same  molecule 
extracted  from  different  crystal  structure  determinations  because 
of  geometrical  distortions  due  to  different  accuracies  in  atomic 
positional  parameters,  which,  however  small,  have  large  con¬ 
sequences  for  the  intramolecular  binding  geometry.  The  effect 
is  less  pronounced  with  intermolecular  energies. 

III.  Statistical  and  Theoretical  Methods 

1.  Atom-Atom  Contact  Distributions.  The  structural 
database  of  Table  1  was  used  to  analyze  atom-atom  contact 
densities.  Distance  distribution  functions  (DDF)  are  defined  as 
follows:8  All  atom— atom  distances  between  atomic  species  i 


TABLE  1:  Compound  Designator  Key:  Cambridge 
Structural  Database  REFCODE,  Explosive  Acronym  (if 
Extant),  and  Common  Chemical  Name 


CATJ1Q 

CIWMEA 

CTMTNA 

DATNBZ 

DIMNAN 

DNBENZ 

DNITBZ 

DNPMTA 

GEMZAZ 

GIMBOT 

HIHHAH 

JEDSUG 

JEHLAJ 

JEXLUT 

JEXMAA 

JEXMEE 

JEXMII 

KEMTIF 

KOFKAR 

KOFKEV 

KOFKIZ 

METNAM 

MTNANL 

MNTDMA 

NACPON 

NITOLU 

OCHTET 

NOHTAZ 

PUBMUU 


SECVOL 

TATNBZ 

TNBENZ 

TNIOAN 

TNOXYL 

TNPHNT 

ZZZFYW 

ZZZMUC 

ZZZQSC 


1 ,3,3,5,7,7-hexanitro- 1 ,5-diazacyclo-octane 
TNAZ,  1,3,3-trinitroazetidine 
RDX,  cyclotrimethylene-trinitramine 

l,3,5-trinitro-2,4-diaminobenzene 
N,N-dimethyl-4-nitroaniline 

1.3- dinitrobenzene 

1 .4- dinitrobenzene 
dinitropenatmethylenetetramine 

1.4- dinitroglycoluril 

2.2.4.4.6.6- hexanitrostilbene 

3- nitrobenzaldehyde 

1 .3- dinitro-l  ,3diazacyclopentan-2-one 

2.5 .7.9- tetranitro-8-oxo-2,5 ,7 ,9-tetra- 
azabicyclo(4.3.0)nonane 

cis-l,3,5,7-tetranitro-l,3,5,7-tetra-azadecalin 
TNBI  l,l,3,3-Tetranitro-4,4-bi-imidazolidine 
TNSD,  1, 3,7, 9-Tetranitro- 1,3,7, 9-tetra- 
azaspiro(4.5  )decane 

trans- 1 ,4,5,8-tetranitro- 1 ,4,5,8-tetra-azadecalin 

2.4.8. 1 0- tetranitro-2,4,8, 1 0-tetra-azaspiro(5  .5)undecane 

1 .3- dinitro- 1 ,3-diazacyclohexane 

1 .3- dinitro- 1 ,3-diazacycloheptane 

1 .5- dinitro-3-nitroso-l  ,3,5-triazacycloheptane 
N,N-dimethylnitramine 

TETRYL,  2,4,6-trinitro-N-methyl-N-nitroaniline 

N,N-dimethyl-3-nitroaniline 

N,N-dimethyl-4-nitroaniline 

4- nitrotoluene 

HMX,  1 ,3,5,7-tetranitro- 1 ,3,5,7-tetra-azacyclo-octane 

2.4.6- trimethyl- 1 ,3,5-trinitrohexahydro,  1 ,3,5-triazine 
HNIW  2, 4, 6,8, 1 0, 1 2-hexanitro-2,4,6,8, 10,12- 

hexa-azatetracyclododecane, 

hexanitrohexa-azaisowurtzitane 

2.2- bis(  1 ,3-dinitrohexahydropyrimidine 

1. 3. 5- triamino-2, 4,6-trinitrobenzene 

1 .3.5- trinitrobenzene 

1.3.5- trinitroaniline 

1.3.5- trinitro-2,4-dimethylbenzene 
1  -ethoxy-2,4, 6-trinitrobenzene 

1.2- dinitrobenzene 

2.4.6- trinitrotoluene 

2.4- dinitrotoluene 


and  j  in  all  crystal  structures  in  the  database  are  calculated  up 
to  a  certain  limit  Rmax.  Let  N^R)'!  be  the  number  of  atom  — 
atom  distances  within  the  Mh  distance  bin  between  Rk  and  /?* 
4-  dR.  An  inspection  of  the  unnormalized  distribution  function 
reveals  a  separation  below  which  no  contacts  are  observed,  R°‘l. 
The  normalization  factor,  Fn,  and  the  DDF  are  then 


FNjj  =  2>*(/0ij/[47r/3((*max)3  -  (/?°'03)] 
DDF  =  gk(R)IJ  =  (1/  F  KiJ)Nk(Rf/dVk 


where  dP*  is  the  volume  of  the  spherical  shell.  Fn  represents 
the  condition  of  uniform  distribution  of  the  observed  contacts 
over  the  available  contact  space.  The  DDF  is  formally  similar 
to  a  radial  distribution  function,  RDF,  but  the  physicochemical 
meaning  is  different,  because  a  RDF  refers  to  a  sampling  over 
the  same  molecular  species  in  a  homogeneous  system,  while 
the  DDF  refers  to  sampling  over  different  molecules  and  over 
many  different  systems,  the  crystal  structures. 

2.  Ab  Initio  Calculations.  Quantum  chemical  calculations 
on  isolated  molecules  with  a  fixed  geometry  as  extracted  from 
the  crystal  structure  determination  were  carried  out  using  the 
GAUSSIAN  package,9  at  the  MP2  6-3 1G**  level  with  the 
“cube”  option  to  produce  an  electron  density  for  PIXEL  (step 
of  0.08  A),  and  the  “pop  =  esp”  option  for  the  calculation  of 
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TABLE  2:  Crystal  Lattice  Energies  (kJ  mor1) 

Eckhardt  and  Gavezzotti 

CSD  refcode 

-£(UNiy 

A//(subl)1’ 

—f  (PIXEL/ 

-fCSRTV 

£(coul+pol/ 

£(disp  y 

CATJIQ 

152 

183 

CIWMEA 

93 

84 

113 

-61 

-92 

CTMTNA 

109 

112 

116 

130 

-87 

-97 

DATNBZ 

154 

140 

116 

-75 

-121 

DIMNAN 

108 

103 

105 

-60 

-118 

DNBENZ 

101 

85 

80 

-40 

-77 

DNITBZ 

105 

96 

89 

-49 

-80 

DNPMTA 

130 

130 

141 

-87 

-107 

GEMZAZ 

148 

155 

-157 

-125 

GIMBOT 

205 

180 

HIHHAH 

101 

110 

88 

-49 

-85 

JEDSUG 

93 

112 

124 

-86 

-81 

JEHLAJ 

136 

181 

JEXLUT 

159 

179 

JEXMAA 

154 

176 

JEXMEE 

147 

173 

JEXM1I 

155 

169 

KEMTIF 

160 

184 

KOFKAR 

101 

101 

117 

-68 

-82 

KOFKEV 

108 

106 

120 

-67 

-92 

KOFKIZ 

117 

131 

METNAM 

59 

70 

72 

70 

-41 

-56 

MTNANL 

134 

134 

118 

150 

-77 

-106 

MNTDMA 

104 

93 

94 

-49 

-105 

N  AC  PON 

100 

NITOLU 

86 

79 

71 

-30 

-73 

NOHTAZ 

127 

141 

OCHTET00 

127 

179 

OCHTET03 

124 

168 

OCHTET12 

146 

162 

157 

180 

-91 

-127 

PUBMUU00 

156 

175 

PUBMUU01 

165 

181 

PUBMUU02 

158 

187 

SECVOL 

152 

188 

TATNBZ 

190 

168 

156 

-113 

-137 

TNBENZ10 

105 

107 

TNBENZ13 

107 

107 

86 

-54 

-94 

TNIOAN 

132 

120 

100 

-62 

-97 

TNOXYL 

132 

130 

103 

-58 

-104 

TNPHNT 

119 

121 

104 

-63 

-100 

ZZZFYW 

101 

87 

88 

-46 

-76 

ZZZMUC 

119 

105 

ZZZQSC 

108 

98 

89 

-46 

-91 

<J  Lattice  energy  with  UNI  atom -atom  force  field.1 1  h  Heats  of  sublimation. c  Total  PIXEL  lattice  energy.  d  Lattice  energy  reported  in  ref  17,  with 
atomic  charges  from  highest-level  quantum  chemical  calculation.  *■  Sum  of  Coulombic  and  polarization  PIXEL  terms,  f  PIXEL  dispersion  energy. 


electrostatic-potential-derived  atomic  point  charges.  Alterna¬ 
tively,  these  were  also  calculated  by  the  so-called  rescaled-EHT 
method,10  based  on  a  Mulliken  population  analysis  on  an 
Extended  Hiickel  calculation  with  modified  valence-orbital 
ionization  potentials.  This  produces  atomic  point  charges  very 
similar  to  those  from  a  Mulliken  population  analysis  on  the  MP2 
wavefunction,  and  can  be  routinely  applied  even  to  very  large 
molecular  systems  because  of  the  virtually  null  computational 
requirements  of  an  EHT  calculation  versus  a  full  molecular 
orbital  MP2  calculation.  These  atomic  point  charges  were  used 
only  for  the  calculation  of  cell  dipoles  in  the  convergence 
correction,  or  for  comparison  with  PIXEL  Coulombic  energies 
which  include  penetration  energies  (see  below). 

3.  Lattice  Energies  and  Molecule- Molecule  Energies. 
Crystal  structures  were  taken  as  is  from  the  X-ray  diffraction 
determination  without  optimization  of  lattice  energies  or 
structure  relaxation,  but  hydrogen  atom  positions  were  re¬ 
calculated  as  usual.10  Force-field  calculations  were  carried  out 
in  the  atom— atom  approximation  with  the  chargeless  UNI  force 
field,10  which  uses  a  simple  A  exp (~BR)  -  CR~ 6  functional 
form  in  the  atom-atom  distance  R.  Point-charge  Coulombic 
terms  were  then  evaluated  separately  using  atomic  point-charge 
parameters  obtained  as  described  above.  For  more  accurate 


calculations,  the  SCDS-PIXEL  approach11  was  applied.  The 
molecular  electron  density  is  first  calculated  by  standard 
quantum-chemical  methods  giving  the  electron  distribution  on 
a  large  number  (10  000—20  000  for  our  molecules)  of  charge 
pixels.  The  Coulombic  energy  is  then  calculated  by  sums  over 
pixel— pixel,  pixel-nucleus,  and  nucleus— nucleus  Coulombic 
terms.  A  local  polarizability  is  then  assigned  to  each  pixel,  the 
electric  field  generated  by  pixels  and  nuclei  in  surrounding 
molecules  is  calculated,  and  the  linear  polarization  energy  is 
evaluated.  An  empirical  damping  function,  using  one  disposable 
parameter,  is  introduced  to  avoid  singularities.  The  overlap 
between  molecular  densities  is  calculated,  and  the  exchange 
repulsion  energy  is  evaluated  as  proportional  to  the  overlap 
integral.  Dispersion  energies  between  two  molecules  A  and  B 
are  calculated  as  a  sum  of  pixel— pixel  terms  in  a  London-type 
expression  involving  the  above-defined  distributed  polarizabili¬ 
ties  and  an  “oscillator  strength”  energy,  Eos'- 

£d.sp.aB  =  (-3/4)  XE£’os^)a(a/[(4^°)2(/f/] 

i, A  »',B 

where  the  damping  function  is^(/?)  =  exp[— (D/R,j  -l)2]  (for 
Rij<D ),  where  D  is  an  adjustable  empirical  parameter.  Standard 
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PIXEL  theory  uses  D  =  3.0  A,  but  a  preliminary  survey  showed 
that  all  lattice  energies  of  crystals  considered  here  were 
underestimated  with  respect  to  the  A//SUb.  Since  the  parameter 
is  fully  adjustable  and  the  effect  is  systematic,  a  correction  of 
this  parameter  to  reduce  the  discrepancy  was  considered 
acceptable.  D  =  2.6  A  is  used  here,  which  increases  the 
calculated  dispersion  energies  by  5—10%  with  respect  to  the 
standard  parametrization. 

Eos  may  be  approximated  by  the  molecular  ionization 
potential  for  small  molecules  or  as  the  energy  of  the  highest 
occupied  molecular  orbital  (HOMO),  since  the  interacting 
electrons  are  peripheral  ones  and  hence  are  roughly  at  the 
HOMO  energy  level.  A  more  refined  approach  considers  each 
pixel  as  a  separate  oscillator,  with  a  formal  ionization  potential 
/„  which  in  turn  is  a  function  of  the  ionization  potential,  7°, 
pertaining  to  the  atom  to  whose  basin  the  pixel  belongs,  and  of 
the  distance  between  the  pixel  and  the  atomic  nucleus,  /?,: 

Eos  =  (Oj)'a 

/,  =  7°  exp( -/?/?,) 

The  parameter  is  a  function  of  the  atom  type.  This  “variable- 
ionization”  form  of  the  theory  amounts,  in  fact,  to  using  different 
dispersion  energy  coefficients  according  to  the  different  kinds 
of  interacting  atomic  basins,  and  yields  more  accurate  lattice 
energies  in  organic  crystals.8 

PIXEL  calculations  were  applied  to  obtain  separate  Coulom- 
bic,  polarization,  dispersion,  and  repulsion  contributions  to  the 
lattice  energies,  except  for  very  large  molecules,  where  the 
calculation  of  electron  densities  is  too  expensive,  or  for  crystal 
structures  with  more  than  one  molecule  in  the  asymmetric  unit, 
due  to  a  technical  difficulty  with  the  evaluation  of  polarization 
contributions.  For  crystal  structures  in  polar  space  groups,  when 
the  angle  between  the  molecular  dipole  direction  and  a  polar 
axis  is  small,  a  convergence  correction  should  be  applied  to 
the  lattice  sums  using  the  van  Eijck — Kroon  method;12  for  the 
very  few  cases  when  they  were  needed,  however,  these 
corrections  were  around  1  —2  kJ  mol  ',  so  they  were  considered 
to  be  within  the  computational  noise  and  were  neglected.  More 
information  on  separate  driving  forces  in  crystal  packing  (the 
“structure  determinant”  packing  analysis)  was  then  obtained 
from  PIXEL  interaction  energies  between  near-neighbor  pairs 
of  molecules  in  the  crystal. 

IV.  Results  and  Discussion 

1.  Density  Distribution  Functions  and  Packing  Efficiency. 

Figure  1  shows  that  in  the  database  of  nitro-derivative  crystal 
structures,  the  oxygen— oxygen  intermolecular  distance  distribu¬ 
tion  function  (DDF)  shows  a  well-developed  peak  at  about 
/?(0**0)  =  3.2  A,  a  peak  that  does  not  appear  when  the  same 
survey  is  conducted  on  a  sample  of  ~~30  000  general  organic 
crystals.  We  initially  interpreted  this  unexpected  fact  as  meaning 
that,  due  to  the  relevant  Coulombic  interactions  present  in  these 
crystals,  particularly  those  due  to  stabilizing  interaction  between 
nitro-nitrogen  sites  and  nitro-oxygen  sites,  the  peripheral  oxygen 
atoms  are  forced  into  a  close  intermolecular  contact  that  is 
counterintuitive  in  terms  of  simple  atom -atom  Coulombic 
interactions.  The  explanation  actually  requires  a  subtler  analysis 
including  penetration  energies,  i.e.,  those  attractive  energies 
which  arise  from  overlap  of  diffuse  electron  clouds  and  their 
proximity  to  positively  charged  nuclei  in  the  intermolecular 
partner. 

A  similar,  but  considerably  smaller,  peak  is  observed  in  the 
DDF  of  0***H  distances.  This  is  not  unexpected  in  view  of  the 


Figure  1.  The  density  distribution  function  of  O**0  separations  in 
the  sample  of  crystals  of  energetic  materials  of  Table  1  (“nitros”  line) 
and  in  a  sample  of  27  461  crystals  of  middle-size,  general  organic 
molecules  containing  only  C,  H,  N,  O,  and  halogen  atoms,  from  the 
Cambridge  Structural  Database. 

opposite  sign  of  the  molecular  potentials  around  hydrogen  atoms 
and  oxygen  atoms  in  nitro  compounds,  and  a  similar  trend  had 
been  previously  noted  in  a  systematic  study  of  nitrobenzenes.13 
These  results  show  the  difficulty  of  defining  intermolecular 
bonds  in  tenns  of  short  atom— atom  intermolecular  distances. 
Thus  the  postulate  of  existence  of  a  C— H***0  crystal  bond  on 
the  basis  of  atom-atom  contacts  alone  requires  explanation  why 
the  same  analysis  does  not  lead  to  the  postulate  of  some  kind 
of  0"0  “bond”  which  would  be  in  conflict  with  elementary 
chemical  reasoning.  Otherwise,  supplementary  reasoning  must 
be  produced  for  the  distinction  between  the  two  cases,  an 
alternative  that  is  also  not  straightforward. 

The  energetic  performance  of  materials  is  often  related  to 
high  crystal  density.  This  property  stems  from  two  unrelated 
factors:  (a)  molecular  mass,  or  the  number  of  protons  and 
neutrons  and  (b)  intermolecular  close  packing,  related  to 
electronic  interactions.  Overall,14  crystal  density  is  determined 
mainly  by  molecular  mass,  an  obvious  property.  For  the  design 
of  materials  for  better  transmission  of  heat  and  mechanical 
energy,  factor  b  above,  control  of  the  relative  efficiency  by 
which  space  is  occupied  by  constituent  molecules  of  identical 
stoichiometry,  is  more  challenging.  A  comparison  between  the 
packing  coefficients15  for  the  present  database  and  for  a  sample 
of  general  organic  crystals  (Figure  2),  in  spite  of  the  large 
difference  in  statistical  significance  due  to  the  large  difference 
in  sample  size,  is  convincing  in  this  respect.  This  observation 
is  undoubtedly  related  to  the  relatively  low  number  of  hydrogen 
atoms  in  these  molecules,  which  permits  a  higher  interaction- 
energy  density  and  parallels  the  previous  observation  of  an 
unusual  distribution  of  short  contacts  between  peripheral  oxygen 
atoms 

2.  Lattice  Energies.  An  extensive  set  of  lattice  energy 
calculations  using  an  empirical  atom— atom  force  field  for 
energetic  materials  with  high  quality  ab  initio  atomic  point 
charges  has  been  carried  out  by  Sorescu,  Rice  and  Thompson 
(SRT).16  We  use  the  lattice  energy  calculated  by  the  universal 
UNI  force  field10  as  a  guideline,  and  Figure  3  plots  the  results 
obtained  by  the  PIXEL  scheme  and  by  the  SRT  scheme,  as  well 
as  the  A//sub.  These  data  are  collected  in  detail  in  Table  2.  The 
best  apparent  agreement  between  experiment  and  calculation 
is  that  of  the  UNI  potentials,  with  SRT  potentials  overestimating 
the  lattice  energies  and  PIXEL  values  showing  an  acceptable 
performance  except  for  a  few  patent  outliers  (TNOXYL, 
TNIOAN,  DATNBZ,  HIM  HAH,  see  also  Table  2)  where  the 
calculated  lattice  energy  is  too  small.  One  should  recall, 
however,  that  the  comparison  is  hampered  on  the  one  hand  by 
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packing  coefficient 


Figure  2.  (a)  Histogram  of  the  percent  distribution  of  packing 
coefficients  in  the  general  sample  (see  caption  to  Figure  1).  The  peak 
value  is  0.70.  (b)  The  distribution  of  the  packing  coefficients  in  the 
sample  of  crystals  of  energetic  materials  of  Table  1  with  average  value 
of  0.77. 


UNI  lattice  energy 

Figure  3.  The  UNI  force-field  lattice  energy  is  the  abscissa:  on  the 
ordinate  are  the  sublimation  enthalpies  (black  dots),  the  PIXEL  lattice 
energies  (squares),  and  the  SRT  lattice  energies  (circles);  see  Table  2 
for  the  definitions.  Units  in  kJ  mol-1. 

the  uncertainties  of  experimental  values,  as  previously  discussed, 
and,  on  the  other  hand,  by  the  neglect  of  intramolecular 
rearrangement  energies  on  going  from  the  crystal  to  the  gas- 
phase  conformation.  Accordingly,  the  overestimation  of  rigid- 
body  lattice  energies  might  be  reasonable,  where  the  correct 
sublimation  enthalpy  results  after  subtraction  of  the  intra¬ 
molecular  energies. 

The  issue  of  accurate  reproduction  of  experimental  thermo¬ 
chemical  data  is  complex,  as  previously  discussed,  but  we  argue 
that  the  added  value  of  the  PIXEL  method  consists  in  the  energy 
partitioning  (Table  2),  which  offers  some  insight  into  the  nature 
of  the  energetic  texture  of  crystalline  materials  not  apparent  in 
atom -atom  schemes.  PIXEL  Coulombic  energies  are  parameter- 


free  and  are  as  accurate  as  the  wavefunction  is,  while  PIXEL 
dispersion  energies  depend,  to  some  extent,  on  the  parametriza- 
tion.  However,  since  Coulombic  energies  are  exact  and  total 
lattice  energies  are  quite  comparable  to  experimental  heats  of 
sublimation,  the  dispersion  contribution  cannot  be  grossly  off 
the  mark.  Therefore,  the  message  in  Table  2  is  unequivocal: 
all  lattice  energies  have  a  substantial  Coulombic-polarization 
(what  is  commonly  called  “Coulombic”)  contribution,  but 
dispersion  energies  are  invariably  equal  to,  or  larger  than,  the 
Coulombic  ones.  This  result,  together  with  the  previous  analysis 
of  oxygen— oxygen  contacts,  issues  a  warning  against  modeling 
schemes  that  rely  too  much  on  Coulombic  energies,  especially 
when  they  are  estimated  by  atomic  point-charge  schemes,  even 
for  the  compounds  considered  here  that  are  commonly  thought 
to  possess  a  strong  permanent  polarization  of  their  electron 
clouds. 

3.  Recognition  Modes.  Significant,  further  information  on 
packing  modes  and  packing  forces  in  organic  crystals  comes 
from  a  calculation  of  interaction  energies  between  pairs  of  close 
neighbor  molecules  in  the  PIXEL  approach,  the  “structure 
determinant”  analysis.17  To  do  this,  one  reference  molecule  is 
chosen  and  crystal  symmetry  is  used  to  generate  those  molecules 
whose  centers  of  mass  are  closest  to  the  center-of-mass  of  the 
reference  molecule.  Then,  each  of  these  coordinated  molecules 
is  designated  by  (i)  the  distance  between  its  center-of-mass  and 
the  reference  molecule,  (ii)  the  label  of  the  symmetry  operation 
connecting  the  two  molecules,  and  (iii)  the  molecule-molecule 
interaction  energy  calculated  by  PIXEL  and  apportioned  ac¬ 
cording  to  its  various  contributions:  Coulombic,  dispersion,  and 
repulsion.  A  set  of  these  geometrical  and  energetic  values  for  a 
molecular  pair  is  referred  to  as  a  “structure  determinant,”  and 
a  set  of  structure  determinants  is  an  unequivocal  and  accurate 
fingerprint  of  a  crystal  structure. 

Data  from  this  procedure  are  more  informative  than  overall 
lattice  energies  because  the  PIXEL  energy  partitioning  quan¬ 
titatively  clarifies  the  relative  importance  of  Coulombic  and 
dispersion  terms,  and  establishes  a  more  immediate  correlation 
between  packing  forces  and  molecular  electronic  structure.  Table 
3  collects  these  lists  for  three  representative  crystals,  and  Figures 
4—6  show  the  corresponding  structural  diagrams. 

In  the  methylnitramine  crystal,  the  two  determinants  along 
the  screw  axis,  A  and  C,  roughly  correspond  to  antiparallel  nitro 
groups,  but  an  interpretation  based  only  on  opposing  dipoles  is 
not  warranted,  given  that  the  Coulombic  and  dispersion 
contributions  are  roughly  equivalent.  The  nitro-methyl  B 
determinant  has  a  predominant  Coulombic  character.  The  next 
three  determinants  complete  the  first  coordination  shell  made 
of  12  nearest-neighbors  as  in  the  close  packing  of  spheroids 
(recall  that  each  S  or  T  interaction  has  two  partners  to  the 
reference  molecule).  They  are  rather  weak  and  sometimes  have 
a  destabilizing  Coulombic  component,  like  the  E  interaction 
which  presents  a  confrontation  of  methyl  groups,  but  are 
stabilized  by  a  counterbalancing  dispersion  term. 

In  the  CTMTNA  crystal,  the  strongest  binding  comes  from 
an  interaction  over  a  centrosymmetric  molecular  pair  with  the 
nitro-oxygens  in  one  partner  pointing  to  the  nitro-nitrogen  in 
the  other.  Here  the  Coulombic  term  is  strongly  stabilizing,  as 
revealed  by  the  PIXEL  calculation,  in  spite  of  the  unavoidable 
proximity  of  oxygen  atoms.  Notice,  however,  that  in  this  case, 
too,  the  predominating  factor  is  dispersion.  It  would  have  been 
quite  difficult  to  properly  describe  this  dominating  interaction 
in  pure  atom-atom  terms,  and  only  the  PIXEL  description  is 
adequate.  The  next  determinant,  B ,  is  a  clear  0"H  Coulombic 
interaction.  However,  it  becomes  increasingly  difficult  to  assign 
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determinant  label, 

R(  c.o.m.),  approximate 

symmetry"  EC(ni\h  £poi  £disP  ETcp  Eiol  description 


METNAM 


A 

3.86 

S 

-17 

-3  -16 

+7 

-29 

nitronitro  antiparallel 

B 

6.13 

T, 

-18 

-4  -9 

+9 

-21 

nitro-methyl,  double 

C 

5.36 

S 

-10 

-3  -10 

+6 

-16 

nitronitro  antiparallel, 
offset 

D 

6.58 

T, 

-3 

-1  -5 

+3 
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nitro-methyl,  single 
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S 

+3 

-1  -10 

+3 
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methyl  methyl  offset 

F 

6.07 
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A 
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B 
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Gv 
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14 

-26 

nitro-mcthylcnc 
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7 
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D 
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oxygen— nitrogen 
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F 
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G- 

-5 

-4  -23  19 
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single  0*"H  contact 
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PE 
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Tc 

+5 
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1 

+2 
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°  Label  as  in  Figures  4-6,  distance  between  centers  of  mass  (com) 
(A);  label  of  the  symmetry  operator:  S  =  twofold  screw,  T  = 
translation,  I  =  inversion  center,  G  =  glide  plane. h  PIXEL  partitioned 
energies  and  total  molecule-molecule  energy,  kJ  mol-'. 


Figure  5.  The  structure  determinants  in  the  cyclotrimethylene- 
trinitramine  (CTMTNA)  crystal  structure.  Top  to  bottom:  A,  B,  and 
C  (same  labels  as  in  Table  3).  Dark  caps,  nitro  groups;  white  caps, 
hydrogen  atoms. 


Figure  4.  The  structure  determinants  in  the  N,N-dimethylnitramine 
(METNAM)  crystal  structure.  Top  to  bottom:  A,  B,  C,  and  D  (same 
labels  as  in  Table  3). 

an  unequivocal  label  to  the  subsequent  determinants,  which 
comprise  broad  and  diffuse  mixtures  of  Coulombic  and  disper¬ 
sion  terms.  The  presence  of  these  nonspecific  structure  deter¬ 


p  phave,  E 

Figure  6.  The  structure  determinants  in  the  a  phase  HMX  (OCHTETOO) 
and  P  phase  HMX  (OCHTET12).  Same  labels  as  in  Table  3.  Dark  caps, 
nitro  groups;  white  caps,  hydrogen  atoms. 


minants  is  what  makes  crystal-structure  analysis,  not  to  mention 
prediction  and  control,  so  difficult  in  terms  of  intermolecular 
atom— atom  bonding. 

The  same  analysis  on  the  A  and  D  molecular  pairs  in  the 
CTMTNA  crystal  has  been  carried  out18  with  a  methodology 
that  in  many  respects  resembles  the  PIXEL  approach,  using  a 
HF/6-311G**  wavefunction.  The  use  of  HF  instead  of  MP2 
leads  to  higher  Coulombic  energies,  as  is  well-known.  On  the 
other  hand,  a  comparison  between  the  results  of  the  two  methods 
at  the  same  wavefunction  level  (HF)  gives  very  close  values 
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point  -charge  coulombic  energy 

Figure  7.  PIXEL  exact  Coulombic  energy  vs  approximate  point-charge 
Coulombic  energy  for  the  crystal  structures  listed  in  Table  2.  kJ  mol**1 
units. 

for  the  Coulombic  +  polarization  term  (—35.2  and  —16.8, 
PIXEL,  -37.2  and  -18.8,  ref  18). 

The  case  with  MMX  is  even  more  complex.  Notice  how  the 
same  molecule  in  two  largely  different  confonnations  can  build 
up  crystal  structure  determinants  of  nearly  equal  intensity  in 
two  polymorphic  crystals:  the  a  form  by  a  clear  confrontation 
of  its  nitro  and  methylene  moieties  with  Coulombic  stabilization 
(ol4  determinant),  and  the  /?  form  by  a  less  clear-cut  “nested” 
contact  (J3A  determinant),  with  an  equal  amount  of  Coulombic 
and  dispersive  character.  In  the  /3  form,  a  near-neighbor  pair 
( (3E)  displays  a  sharp  confrontation  of  opposing  nitro  groups 
where  the  distance  prevents  an  overlap  of  the  molecular  electron 
densities  and  no  penetration  energy  sets  in.  The  resulting 
Coulombic  interaction  energy  is  mildly  destabilizing,  but  a  minor 
amount  of  dispersion  from  the  overall  molecular  bodies  is  more 
than  enough  to  wipe  out  the  repulsion.  Thus,  the  PIXEL  analysis 
reveals  the  presence  of  nonspecific  stabilizing  or  even  “silent”, 
slightly  destabilizing  partners  whose  role  in  crystal  packing  is 
difficult  to  assess  but  is  by  no  means  negligible,  although  they 
are  usually  neglected  in  standard  crystal  packing  analyses  based 
on  short  interatomic  distances  only. 

4.  Point-Charge  versus  PIXEL  Coulombic  Energies.  A 
comparison  between  PIXEL  Coulombic  lattice  energies  and 
point-charge  Coulombic  energies  is  instructive.  Figure  7  shows 
the  result:  the  good  news  is  that  there  is  a  rather  consistent 
proportionality  ratio  by  which  the  point-charge  energy  is  about 
60%  of  the  real  Coulombic  energy.  This  could  be  useful  in  the 
rescaling  of  empirical,  force-field  parametrizations.  There  are, 
however,  some  patent  outliers  where  the  point-charge  model  is 
severely  off  the  mark.  The  case  of  triamino-trinitrobenzene,  with 
its  extensive  hydrogen-bonding  network,  is  exemplary  in  this 
respect. 

V.  Conclusions 

The  quantitative  analysis  of  structures  and  interaction  energies 
in  crystals  of  energetic  compounds  is  in  many  respects  revealing. 
With  their  high  number  of  nitro  groups,  these  molecules  have 
a  low  relative  abundance  of  hydrogen  atoms  and  in  their  crystals 
they  must  make  do  with  an  unusually  high  number  of  contacts 
between  peripheral  oxygen  atoms.  These  are  indeed  numerous 
and  anomalously  short,  as  clearly  appears  in  a  quantitative  study 
based  on  the  newly  introduced  concept  of  Density  Distribution 
Functions  (DDF’s).  Packing  coefficient  values  are  unusually 
high,  most  likely  as  a  consequence  of  a  higher  interaction  energy 
density.  Calculated  lattice  energies  compare  rather  favorably 
with  experimental  sublimation  enthalpies  in  all  of  the  compu¬ 
tational  schemes  employed  here.  The  comparison  is,  however, 
much  less  compelling  than  is  usually  thought  in  molecular 


simulation  and  force-field  calibrations  because  of  the  large 
uncertainties  in  experimental  values  and  because  of  the  many 
approximations,  especially  the  rigid-molecule  approximation, 
which  is  particularly  restrictive  for  many  of  the  compounds 
examined  here.  The  PIXEL  analysis  has  an  average  performance 
in  regard  to  total  lattice  energies  but  is  superior  in  its  interpreta¬ 
tive  power  because  it  is  based  on  complete  electron  distributions 
rather  than  on  a  few,  localized,  interaction  centers.  The  PIXEL 
analysis  reveals  that  (i)  short  oxygen-oxygen  separations  do 
not  imply  repulsive  or  destabilizing  interactions,  a  probable 
conclusion  arising  from  a  superficial  analysis  in  terms  of  atom- 
atom  contacts,  and  (ii)  dispersive  contributions  to  the  interaction 
energies  are  invariably  higher,  contrary  to  common  belief. 

The  studies  presented  here  use  molecular  simulations  based 
on  a  perfect  model,  an  infinite,  homogeneous  crystal  without 
defects.  Whether  this  model  is  adequate  or  not  depends  on  which 
crystal  properties  are  studied.  We  suggest  that  our  results  are 
significant  concerning  the  general  conclusions  regarding  struc¬ 
tural  features  of  these  crystals:  atom -atom  distance  distribu¬ 
tions  and  packing  coefficients,  and  the  general  nature  of  the 
crystal  packing  forces,  especially  their  subdivision  into  Cou¬ 
lombic  and  dispersive  components. 

In  any  case,  simulation  of  perfect  crystals  is  expected  to 
provide  a  paradigm  by  which  the  effect  of  imperfections  can 
be  discussed.  Moreover,  molecule— molecule  interaction  energies 
do  not  depend  on  a  global  model  of  the  bulk  crystal  and  are 
therefore  certainly  informative  and  reliable.  A  proper  consid¬ 
eration  of  the  role  of  defects  will  certainly  be  needed  when 
tackling  a  molecular  simulation  for  a  more  detailed  investigation 
of  subtle  factors  which  determine  the  properties  related  to 
reactivity,  decomposition  and  mass  or  heat  transport  of  these 
materials.  This  is  left  for  future  work. 
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